% ***********************************************************************************************
% m_no_agglomeration.m 
% compute equilibrium when we shut down agglomeration in productivity and amenities 
% ***********************************************************************************************
clear all; 
close all;
format shortG;
% ==================== LOAD FUNDAMENTAL PARAMETERS ========================
load parameters.mat;                            
load calib_productivity_amenity.mat;    
load calib_setup.mat;                    
load commuting_cost.mat;  
load calib_auxiliary.mat;

load fundamentals_1950.mat;

load calib_qadjusted.mat; 
load calib_hscale.mat; 

Rdata=Rdata_all; Ldata=Ldata_all;  
totalpop=sum(Rdata,1);

N=size(ID,1);
NT=6; 
% ==================== COMPUTE COMMUTING PATTERN 1945 =====================
Lcom=zeros(N,N,NT+1); 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-5,'TolX',1e-5);
R45=Rdata(:,1); L45=Ldata(:,1); 
K45=Kappa_mat(:,:,1).^(-rho/sigma); 
Y0=ones(N,1);     
comsyst = @(y)fcomeval2(y,K45,R45,L45,rho,sigma);
[y_fsolve, yfval_fsolve]=fsolve(comsyst, Y0, options0); Y45=abs(y_fsolve); Y45=Y45./geomean(Y45);  
ytemp=K45.*(repmat(Y45,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Lcom45=ytemp.*repmat(R45',N,1); 
Lcom(:,:,1)=Lcom45;

% ==================== EXOGENOUS FUNDAMENTALS (50-75) =====================
a_mat=[a50, aa];   a_mat=a_mat./repmat(geomean(a_mat),N,1);
b_mat=[b50, bb];   b_mat=b_mat./repmat(geomean(b_mat),N,1);

% ================= COMPUTE OPTION VALUES OF FUNDAMENTALS =================
Xmat_noagg=zeros(N,NT); Ymat_noagg=zeros(N,NT); 
Xmat_noagg(:,NT)=b_mat(:,NT); Ymat_noagg(:,NT)=a_mat(:,NT); % last period (1975)
for t=1:NT-1
    tt = NT-t; 
    theta=thetavec(tt);
    Xtemp=Xmat_noagg(:,tt+1).^(rho*(1-theta)); Xtemp=b_mat(:,tt).*Xtemp;
    Ytemp=Ymat_noagg(:,tt+1).^(rho*(1-theta)); Ytemp=a_mat(:,tt).*Ytemp;
    Xmat_noagg(:,tt)=Xtemp;
    Ymat_noagg(:,tt)=Ytemp;
end

%% ================= COMPUTE COUNTERFACTUAL FROM 1950 ======================
Q_i=ones(N,NT); 
L_i=Ldata(:,2:end); R_i=Rdata(:,2:end); % initial guess 

Q_n=zeros(N,NT); L_n=zeros(N,NT); R_n=zeros(N,NT); 
H_n=zeros(N,NT+1); H_n(:,1)=hsupply45; 
H_demand=zeros(N,NT); 

maxit=1e+6; tol=1e-6; update=0.05; it=0; 

while it<maxit  
   
LQmat=zeros(N,NT); 
LQmat(:,NT)=Q_i(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    LQmat(:,tt)=Q_i(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta)));
end 
 
for t=1:NT
    if t==1; theta=theta50; Rpre_i=R45; Lpre_i=L45; 
    else     theta=thetavec(t-1); Rpre_i=R_n(:,t-1); Lpre_i=L_n(:,t-1); end 
    
    lambda_L=Kmat(:,:,t).*((repmat(LQmat(:,t)',N,1).^(-mu)).*(repmat(Xmat_noagg(:,t)',N,1))).^(rho/sigma); 
    lambda_L=lambda_L./repmat(sum(lambda_L,2),1,N); 
    lambda_R=Kmat(:,:,t).*((repmat(LQmat(:,t),1,N).^(-gammaH/gammaL)).*(repmat(Ymat_noagg(:,t),1,N).^(1/gammaL))).^(rho/sigma);
    lambda_R=lambda_R./repmat(sum(lambda_R,1),N,1);
 
    Rpost_i=R_i(:,t); Lpost_i=L_i(:,t);
    it1=0;
    while it1<maxit 
 
    Rpost_n=(1-theta).*Rpre_i+lambda_L'*(Lpost_i-(1-theta).*Lpre_i);
    Lpost_n=(1-theta).*Lpre_i+lambda_R*(Rpost_i-(1-theta).*Rpre_i);  

    if max(abs(Rpost_n-Rpost_i))<tol && max(abs(Lpost_n-Lpost_i))<tol; break;
    else
       Rpost_i=update.*Rpost_n+(1-update).*Rpost_i; 
       Lpost_i=update.*Lpost_n+(1-update).*Lpost_i;
       it1=it1+1;  
    end
    end
    R_n(:,t)=Rpost_n; L_n(:,t)=Lpost_n; 

    % commuting patterns 
    Lcompost=(1-theta).*Lcom(:,:,t)+lambda_R.*repmat((Rpost_n-(1-theta).*Rpre_i)',N,1); 
    Lcom(:,:,t+1)=Lcompost;
    % wage (workplace) and average income (residential place) 
    w_n=(Q_i(:,t).^(-gammaH/gammaL)).*(a_mat(:,t).^(1/gammaL));
    v_n=(Lcompost./repmat(Rpost_n',N,1))'*w_n; 
    % aggregate demand and supply (stock) of floor space
    h_n=(1-psi).*H_n(:,t)+hbar(t)*(Q_i(:,t).^eta).*Area;
    H_n(:,t+1)=h_n; 
    floor_demand=mu.*v_n.*Rpost_n + (gammaH/gammaL).*w_n.*Lpost_n;
    H_demand(:,t)=floor_demand; 
    % floor space market clearing condition
    Q_n(:,t)=floor_demand./h_n ;
end 

dd_Q=Q_n-Q_i; 

if max(max(abs(dd_Q))) < tol 
        disp(['Converged floor prices']); break;
else
   Q_i=update.*Q_n+(1-update).*Q_i; it=it+1;  
   if rem(it,10)==1
      disp(['Largest Error is:',' ',num2str(it), ' ',num2str(max(abs(dd_Q)))]);
   end
end

end

max(H_n./Area)
median(H_n./Area)

if max(abs(Q_n.*H_n(:,2:end)-H_demand))<1e-6; disp(['PASS floor space market clearing conditions in 1950-75']);
else; disp(['WARNING']); end 

% *******************************************
%       SAVE NO AGGLOMERATION RESULTS
% *******************************************
R_noagg=R_n; L_noagg=L_n; Q_noagg=Q_n; H_noagg=H_n; Lcom_noagg=Lcom; 

result=array2table([R_noagg,Rdata,L_noagg,Ldata, Q_noagg, Qadj, hstock, H_noagg, cbddist, Area, ID, lat, lon],'VariableNames',{'pop50_noagg',...
    'pop55_noagg','pop60_noagg','pop65_noagg','pop70_noagg','pop75_noagg',...
    'pop45','pop50','pop55','pop60','pop65','pop70','pop75',...
    'emp50_noagg','emp55_noagg','emp60_noagg','emp65_noagg','emp70_noagg','emp75_noagg',...
    'emp45','emp50','emp55','emp60','emp65','emp70','emp75',...
    'q50_noagg','q55_noagg','q60_noagg','q65_noagg','q70_noagg','q75_noagg',... 
    'q50','q55','q60','q65','q70','q75',... 
    'h50','h55','h60','h65','h70','h75',... 
    'h45','h50_noagg','h55_noagg','h60_noagg','h65_noagg','h70_noagg','h75_noagg',...
    'cbddist', 'area', 'geocode', 'lat', 'lon'}); 
writetable(result,'../../output/quant/output_no_agglomeration.csv');

save('noagglomeration.mat','R_noagg','L_noagg','Q_noagg','H_noagg','Lcom45','Lcom_noagg');
% ********************************************
%   FUNCTIONS 
% ********************************************
function[comeval] = fcomeval2(Y,K,DR,DL,rho,sigma)

N = size(DR,1); 
Y=abs(Y);  
Y=Y./geomean(Y);
ytemp=K.*(repmat(Y,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
comeval = DL - sum(ytemp.*repmat(DR',N,1),2); 
comeval = abs(comeval); 

end

